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In equilibrium molecular dynamics, Einstein relation can be used to calculate the 
»- Eh ' thermal conductivity. This method is equivalent to Green-Kubo relation and it does 

not require a derivation of an analytical form for the heat current. However, it is not 
commonly used as Green-Kubo relationship. Its wide use is hindered by the lack of a 
proper definition for integrated heat current (energy moment) under periodic bound- 
ary conditions. In this paper, we developed an appropriate definition for integrated 
heat current to calculate thermal conductivity of solids under periodic conditions. 
We applied this method to solid argon and silicon based systems; compared and 
contrasted with the Green-Kubo approach. 
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I. INTRODUCTION 



Accurate determination of thermal conductivity of materials is of fundamental importance 
especially for technological applications where the control of thermal transport is critical. 
For instance, microelectronics, where the device failures are mostly due to overheating, re- 
quire a strict thermal management. The development of efficient thermoelectrics require 
materials and materials systems with low thermal conductivity. A similar requirement arise 
in the development of thermal barrier coatings. Experimental evaluation of thermal con- 
ductivity is not as straightforward as electrical conductivity. Measurements generally suffer 
from heat losses, non-uniform heating and errors introduced due to approximations that 
account for sample size and structure^. As the dimensions of the samples get smaller, these 
effects become more pronounced, increase the inaccuracy of the measurements. Due to 
their wide use in electronics, thermal transport in silicon and silicon based alloys has been 
studied extensively. Within the last decade, silicon based materials have been considered 
for thermoelectric and thermal barrier applications where a small thermal conductivity is 
desired 2 . The studies in this area generally focus on the influence of defects, impurities and 
nano-structure (superlattices^, nanowires^, quantum dots^, etc.) on thermal transport. 

As an alternative to experiments, computational approaches offer controllable virtual 
measurements to determine thermal conductivity, with the added potential of probing the 
processes and mechanisms that influence the thermal transport. The two most frequently 
used computational methods that include atomistic level details in studying thermal trans- 
port are based on Boltzmann transport equation (BTE) and molecular dynamics (MD) 
simulations. The first one, BTE, is a semiclassical approach. The atomistic level granularity 
in BTE, is provided by the determination of phonon frequencies and group velocities from 
the lattice dynamics of an N-body system. However, the detailed physics of the heat trans- 
fer is incorporated into the formulation through phonon relaxation times^ - — . This requires 
either experimental results or a detailed analysis of the fundamental scattering processes to 
provide reliable theoretical estimates of relaxation times. When such data is not available 
due to the complexity of the system, the relaxation times can be probed with MD. Details of 
such calculations are described in earlier studies where MD was used to obtain the phonon 
relaxation times in fee crystals^ and silicon^ii using the decay of phonon mode amplitudes 
or energies. 
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In contrast, the full atomistic detail of the physics of heat transport is inherently present 
in MD through the definition of interactions between the atoms that make up the system. 
The evolution of the state of the system is given by the equations of motion based on inter- 
action potentials that are usually anharmonic. Hence, MD is a powerful tool, as it includes 
normal, Umklapp and higher order phonon-phonon scattering events through the dynamics 
of the thermodynamic state of the system. Of course due to the classical nature of MD, 
vibrational modes are equally excited (i.e. have the same energy /cgT)," and the phonon 
scattering processes do not include electronic contributions. However, some approaches" 
have also implemented quantum corrections. There are three MD based approaches for cal- 
culating thermal conductivity: the steady state non-equilibrium (direct) method, synthetic 
hamiltonian non-equilibrium method and the equilibrium (Green-Kubo) method. The direct 
method mimics experiment by imposing a temperature gradient on the system and deter- 
mines the thermal conductivity from the phenomenological Fourier lawJ^-. Temperature 
gradients can be created by scaling the particle velocities, changing the ensemble averages 
or switching particle velocities between the hot and cold ends. There can be some complica- 
tions with this method due to the nonlinear response of a system associated with very large 
temperature gradients. The second technique, proposed by Evans^ applies a homogenous 
synthetic field to create a heat flux, as opposed to applying a temperature gradient. Finally, 
Green and Kubo22^ showed that through the fluctuation-dissipation theorem, transport 
properties i.e. the diffusion coefficient (D), shear viscosity (r)^) and thermal conductivity 
(k) can be calculated by integration of autocorrelation functions of particle velocity (v l ), 
pressure tensor (P^u) and heat current (J). 



where and v are the vector components, V is the system volume, is the Boltzmann 
constant, T is temperature and N is the total number of particles. The same transport 
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coefficients may also been obtained via 'Einstein' relation; which are related to doubly 
integrated form of these expressions. 
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In the Eqs. HI [5]and|6l r l , pi^ and R are the position and momentum of i th particle and 
integrated heat current respectively. The integrated heat current will be referred to as the 
energy moment in this paper and is defined for any N-body system as: 



where is the position and is the energy content of the i particle. 

Detailed derivations on the time correlations functions and associated Einstein relations 
for transport properties can be found in several sources^ - — . In the case of thermal conduc- 
tivity, the heat current autocorrelation function (HCACF) approach has been employed for 
a number of materials such as bulk diamond 14 , bulk silicon^, metal organic frameworks^, 
nanotubes^, nanoribbons^i and nanoparticles 32 . However, unlike the wide use of the Ein- 
stein relation in determining mass transport/diffusion, the use of the Einstein relation for 
thermal conductivity is rather rare. This situation has arisen from a lack of proper imple- 
mentation in the calculation of integrated heat current /energy moment in simulations. 

In the application of the Einstein relation based thermal conductivity calculation in MD, 
discrepancies arise regarding the calculation of the energy moment and energy current. The 
techniques used for calculating R give an expression which produces a bounded behavior and 
results in zero thermal conductivity in non-diffusive solid systems, while J yields a reliable 
result. Similar problems have been reported for the calculation of shear viscosity from the 
Einstein relation using MD with periodic boundary conditions'^ 4 -. Here, we propose an 
expression for R that gives a proper thermal conductivity and can be reduced numerically 
to give the exact value of J. Further, for solid systems, the resulting energy moment 
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expression can be generalized to be independent of the potential model. In this report, 
we formulate the new method and demonstrate its use through applications to a two-body 
potential (Lennard- Jones) and then an N-body potential (Tersoff)^. 

The paper is organized as follows: In the next section we present the background and 
reformulated form of the energy moment based method of calculating thermal conductivity. 
In section III, we present the application of the method, to Argon, Silicon, and Porous 
Silicon and silicon Germanium nanodots in Silicon matrix, providing a detailed comparison 
of energy moment based method with HCACF based methods. Where we observe that 
HCACF of heterogeneous systems (Si-Ge) or systems with extended defects (Si with pores) 
are showing large fluctuations that one needs to filter out or large variance in resulting 
thermal conductivity values. In contrast the Einstein relation for all systems has a much 
better behavior. In section IV, we summarize our results and present our conclusions. 



II. BACKGROUND AND THEORY 

The fundamental variable for calculating thermal conductivity from Eq. [6] is the energy 
moment, R is defined above in Eq. [7J The energy content of particle i, e;, is defined as: 

€i = ^rriivf + Ui (8) 

Here, the first and the second terms are the total kinetic and potential energy of the 
particle i. The velocity and the mass of the same particle are represented by v i and m^. The 
heat current, J, is the time derivative of R. 

As the general definition of R and J are given, we will describe the method of determining 
thermal conductivity with two-body potential and generalize to n-body potential. For two- 
body interactions, the total potential energy content of a particle is defined as Ui = | Y^j^i u ij 
where is the potential between particle i and j. Then, the corresponding heat current 
can be calculated through the following expression: 

- 1 

J = X>^+2 E r dftj-Vi) (10) 

i=l i,j i^j 



5 



where is defined as the nearest image distance and the interatomic force on % due 
to j. 

A straightforward implementation of Eq. [6] through Eq. [8] we call formulation 1, by 
considering a bulk non-convective solid in a periodic box, as we will show produces no 
thermal conductivity. Moreover, if in Eq. [10] is the nearest image distance, then Eq. [7] 
cannot be reduced to Eq. [101 The issue of defining R in periodic systems has been pointed 
out by Donadio and Galli^, who attributed the inconsistency to the ill- definition of particle 
positions in a periodic system. However, we believe this is a consequence of the omission of 
what we call cross-boundary interactions, interactions or transfer of energy between the N 
atoms in the simulation region and the image atoms across the periodic boundary. These 
interactions are vital to accurately defining the energy dispersion of the system and must 
be included into the definition of R. With the inclusion of the cross-boundary interactions, 
energy moment can be separated into potential (Rp) and kinetic (Rk) contributions. The 
combination of R p + Rk will be called formulation 2 for the rest of the report. 

The new form for the potential energy portion of R, Rp, is 



where Uij is the pair potential, i represents atoms in the box and j represents the (real or 
image) neighbors of i. To avoid double counting, j > % constraint is imposed. Though the 
form is very similar to Eq. [3, if an atom in the simulation box interacts with an image atom, 
the actual atomic position of the image is used, instead of contributing the image interaction 
to its corresponding particle in the box. Although a simple equipartition of potential energy 
is the reasonable choice in the case of two-body interactions, n-body interactions may require 
other treatments and a rule of energy distribution may not preferable to another. In these 
circumstances, one needs to derive a separate expression for the heat current from R for each 
choice of potential splitting. However, this derivation is avoided when using R p directly. 

The kinetic portion of the energy moment requires a similar approach. In short, the 
transfer of kinetic energy from atoms in the box to image atoms must be included. This 
transfer is mediated by power term, / • v. The kinetic portion, R^., can then be written as: 




(11) 
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where represents the force on atom % due to j. This equation represents the transfer of 
kinetic energy from atom i to all particles that interact with it. Again, if the neighbor is an 
image atom, the image position is used. 

This discussion has been carried out for the two body interaction, but it can also be 
formatted for n-body potential. In fact, the kinetic part can be put in a form that does not 
depend explicitly on the potential form by simply rearranging the summation as: 



where I is a summation over particles (real or image) that make up unique n-body interaction 
groups. Though similar to the previous expression, Eq. [T3]is well defined for any interatomic 
potential as fi is simply the force on particle I due to all unique interaction groups that 
include 

The computations can be restricted to the kinetic portions of the energy moment when 
dealing with the non-convective systems. In solids, there is very little convection, and 
to see diffusion in a simulation is highly improbable. This, along with the fact that the 
potential contribution to the atomic energy is bound, leads to the conclusion that R p does 
not contribute to the thermal conductivity in perfect solid systems. Then, Rk must be 
main component driving the increase in the Einstein relation. On the other hand, for highly 
diffusive systems such as liquids, R p is expected to contribute significantly. 

III. RESULTS AND DISCUSSIONS 

A. Model Systems 

In order to further illustrate the formulations given in the theory, we studied archetyp- 
ical two-body (argon) and n-body (silicon) systems. Moreover, we discuss the behavior of 
Einstein relation when impurities (germanium) and pores are present in silicon. 

1. Two-Body Potential 

We calculate the lattice thermal conductivity of face centered cubic argon. We have 
chosen the same system size (256 atom) and equilibration procedure as in McGaughey and 
Kaviany^ in order to calculate HCACF. After equilibration, the system is allowed to run 




(13) 
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in NVE ensemble at 50 K for 2 ns with a time step of 0.1 femtosecond. We are aware that 
this time step is small but we chose it to be consistent with the previous literature^ and 
to enable us to calculate more accurate numerical derivatives of (R). We performed five 
independent simulations to average HCACF's. The normalized HCACF is plotted in FigfTJ 
which is similar to what was obtained in McGaughey and Kaviany's work. 

Is 




2 3 
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FIG. 1. (Color online) Normalized HCACF obtained for 256 atom FCC argon in this work is 
represented by (black) solid line. The same function, calculated by McGaughey and Kaviany 



(Ref. 



371 ) is represented by (red) dotted line. 



In Fig. [2] we have plotted the time evolution of R. It is clear that there is a quantitative 
difference between formulation 1 and formulation 2. The former, as expected, fluctuates 
around a constant value. Whereas the latter displays an increasing tendency. 

In principle, the analytical and numerical derivatives of R should produce the same heat 
current («7). Fig. [3] compares the behavior of the heat current obtained from the analytical 
and numerical derivative of R. The inconsistency between the formulation 1 and analytical 
J is obvious. On the other hand, the numerical derivatives of formulation 2 perfectly 
matches the analytical form. 

For the calculation of thermal conductivity, the system size is enlarged to 4000 atoms 
and the total simulation time is extended to 20 ns. Again, five different initial conditions 
are evaluated and the results are averaged. Ensemble average for the autocorrelation is 
obtained by selecting 6 millions time origins 3 fs apart from each other. In calculation of k 
from HCACF, the data was fitted to a double exponential, Eq. [HJ and the infinite integration 
was carried out over this function, an approach adopted by many previous studies^2§»22»21 



S 



-160 
-170 
^-180 

> 

A -190 
Pi 



-200 
-210 



h 



■ AM V v 



v , L_ 

500 



500 1000 1500 

Time (picosecond) 



2000 



FIG. 2. (Color online) Comparison of energy moments (R) for formulation 1, (black) dashed line, 
and formulation 2, (green) solid line. 
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FIG. 3. (Color online) Heat current obtained from analytical, (black) solid line, and numerical, 
(green) dashed line, derivatives of formulation 1. (Grey) dotted line represents numerical derivative 
of formulation 2. 



f(t) = ai e-'/ T1 + a 2 e-'/ T2 (14) 

In this equation oi, a,2, T\ and r 2 are parameters to be fitted and t is the correlation time. 
Both the autocorrelation function and corresponding double exponential function are plotted 
in Fig. HI The thermal conductivity is calculated as 0.403±0.01 W/mK which is very close 
to the value (0.417 W/mK) given in Ref . S 

The Einstein's relations corresponding to formulation 1 is presented in Fig. [5j The 
formulation in Eq. [6] shows that in order to obtain a thermal conductivity, there should 
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FIG. 4. (Color online) Heat current autocorrelation function, (black) solid line, obtained for 4000 
atom FCC argon at 50 K and the corresponding two exponential fit, (green) dashed line. The inset 
shows the same functions in the same units but for shorter correlation time i.e. 10 ps. 



be a finite positive slope of the given function. In the case of formulation 1, two different 
behavior can be identified. First, there is a fast increase in the function and then a plateau 
is reached with no further change thereafter, see the inset of Fig. [5j Due to this behavior, an 
ambiguity arises in the calculation of thermal conductivity. As time goes to infinity the total 
slope goes to zero giving a zero conductivity. Different numbers of time origins from 100 
thousands to 6 millions have been tried to calculate ensemble average for Einstein's relation 
however the behavior of the curve did not change. 
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FIG. 5. Einstein relation using formulation 1 obtained for 4000 atom FCC argon crystal at 50 K. 
The inlet shows the same function in the same units but for shorter time i.e. 40 ps. 
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The second model, on the other hand, has a finite positive slope as seen in Fig. [6j This 
result has been obtained by using 200 thousands time origins with 90 fs apart. In this 
figure, we also present the corresponding linear fit to Einstein relation and the upper-lower 
bounds (corresponding to five different initial conditions). The thermal conductivity from 
the linear fit is 0.416±0.05 W/mK which is very close to the value (0.403 W/mK) obtained 
by HCACF. The variance in the calculated values can be further reduced by using more 
time origins in a longer total run time. 

30000 1 1 1 . 1 . 1 




Time (picosecond) 

FIG. 6. (Color online) Einstein relation using formulation 2, (black) solid line, obtained for 4000 
atom FCC argon crystal at 50 K and the corresponding error limits, (green) dashed lines. (Red) 
dotted line represents the linear fit Einstein relation. 

In the previous section we pointed out that bmR is composed of a kinetic (Rk) and a 
potential {R p ) portion and further stated that R p should not contribute to the thermal 
conductivity of the non-convective solid system. For a diagonal element, Einstein relation 
can be separated and written as 

([R(t) - R(0)] 2 ) = 

([R k (t) - R k (0)} 2 ) + ([R p (t) - R p (0)} 2 ) + 

2([R p (t) - R p (0)][R k (t) - R k (0)]) (15) 

We calculated the first and second terms on the right hand side of the Eq. [15] and plotted 
in Fig. [7J It can be seen that Einstein relation for R k is equal to the Einstein relation for 
total -R. The second term on the right hand side of the Eq. [15] is just a flat line and do 
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not add to the total. Accordingly, the coupled term in Eq. [15j on the average, does not 
contribute to the thermal conductivity. 




Time (picoseconds) 



FIG. 7. (Color online) Einstein relation for kinetic (Rk), (black) solid line, and potential (R p ), 
(red) dashed line, contributions to total R obtained for 4000 atom FCC argon crystal at 50 K. 



2. N-Body Potential 

We selected silicon system represented by the Tersoff potential for many body simulations. 
We applied the same procedure to calculate Einstein relation for 10x10x10 supercell of Si, 
see Fig. [HI By also simulating 6x6x6 and 8x8x8 we have established that the calculated 
conductivity is well converged within 5%. The calculated thermal conductivity of silicon 
is 160.5±10.0 W/mK which is approximately 10% different than the experimental value of 
142 - 148 W/mK — ^ at 300 K. Again, the Einstein relation for Rk is almost identical to 
the one for total R. This confirms that the thermal conductivity does not depend on how 
the potential energy is split between interacting atoms and only the phonons are relevant 
to calculations. Ladd et a/.— draw a similar conclusion that the heat current should be 
independent of the localization of the potential energy in a system at or near equilibrium. 
In their study, Schelling et a/.— re-derived the heat current for different ways of energy 
distribution for a 3-body potential and calculated the thermal conductivity of silicon by 
using HCACF. Their results also did not differ for the reasons stated here. However, splitting 
may be important for the system where convective motion is dominant. 
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FIG. 8. (Color online) Kinetic (Rk), (black) dashed line, and potential (Rp), (red) dotted line, 
portions of total energy moment, (grey) solid line. The data is obtained for 8000-atom silicon 
crystal at 300 K. 

B. Porosity and Ge-Clusters in Silicon 

Heterogeneous systems contain a mixture of elements having different masses, elastic con- 
stants and lattice parameters. Due to this fact, HCACFs show high frequency fluctuations 
through its decay in heterogeneous materials. These high frequency fluctuations pose com- 
putational problems when thermal conductivity is calculated. In order to reveal the true 
decay time, HCACF must be filtered. This is generally done by transforming the function to 
frequency space and eliminating the high frequencies. Sometimes it is necessary to remove 
so many data that the back-transformed function becomes jagged and a good fit to a decay- 
ing function cannot be obtained. This problem is encountered very frequently because the 
heterogeneous materials are studied more often for their technological importance compared 
to their pure counterparts. Being heterogeneous systems, silicon-germanium alloys have the 
same problem when thermal conductivity is evaluated via HCACF. In order to point out the 
benefit of using Einstein relation in heterogeneous systems, we have calculated the thermal 
conductivity of silicon matrix which is embedded with Ge nano-clusters. Each cluster made 
up of 22 germanium atoms surrounded with a regular octahedron of silicon (approximately 
1.65 nm in diameter), see Fig fTDT a). The surfaces of the octahedron consist of {111} family 
of planes. The clusters are placed in a 12x12x12 supercell of silicon and they form a peri- 
odic simple cubic arrangement having a side of 3.311 nm. In Fig[9]^a) we plotted the noisy 
HCACF obtained for the given nano-composite. Without filtering the high frequencies it is 
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very difficult to fit this function to an exponential decay. On the other hand, in Fig. Mjo) it 
can be seen that the Einstein relation behaves much better. The reason of this can be found 
in the integrated nature of Einstein relation where high frequency fluctuations in HCACF 
cancel out. The thermal conductivity is calculated as 21.74 W/mK for this system. This is 
an 80% decrease over the value for pure silicon with only 1.3 at.% Ge impurity. The clusters 
introduce internal surfaces where the opposite sides have different elastic constants, atomic 
masses and lattice parameters. Consequently, processes such as boundary and impurity 
scattering are heavily contributing to the thermal resistance of the lattice. 
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FIG. 9. (Color online) (a) HCACF and (b) Einstein relation for silicon-germanium nano-composite 
having germanium clusters of size 1.65 nm. The clusters form a simple cube having a side length 
of 3.311 nm in a 12x12x12 supercell silicon matrix. 

The minimum thermal conductivity in bulk silicon was estimated as 1.0-1.3 W/mKM-^ 
when the structure is amorphous. These studies suggest that for amorphous silicon the 
phonon mean free path is on the order of average nearest neighbor distance. Later experi- 
mental studies claimed that even lower thermal conductivities than the ones in amorphous 
systems are possible in doped crystalline silicon with nanopores^ 1 ^ 1 ^. This modification, to 
lower extend (i.e. an order of magnitude decrease from bulk thermal conductivity) was also 
observed in microporous silicon by Song and Client. This result is not surprising because in 
addition to decrease in phonon mean free paths and relaxation times due to internal surface 
scattering, there are less conduction channels due to missing silicon atoms. We have tested 
these ideas on porosity with the corrected R and Einstein relation. The porous silicon 
structures and the corresponding thermal conductivities are given in Fig. [TU] and Table |T] 
respectively. In Table HJ a* represents the same octahedron in Fig. [10] without germanium 
atoms inside. Porosity is defined as the ratio of number of missing atoms to number of 
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atoms in a perfect system. We have also presented the ratio of 2 bonded atoms to 4 bonded 
atoms (r2) and the ratio of 3 bonded atoms to 4 bonded atoms (r3) as an indication of the 
internal surface area. The tubular pores in Fig. [TuT b).(c) and (f) have {110} surface planes, 
Fig. [TuT d) has both {100} and {110} surface planes and Fig. [TOT e) has {100} surface planes. 
In general, as the porosity and internal surface area increase, the thermal conductivity de- 
creases as seen in Table |T] in accordance with the previous findings from other theoretical 
methods iL2£. However, this trend breaks down for structures d, e and f. Although f has 
higher porosity and surface area than d and e, it has higher thermal conductivity. Structures 
d and e have atoms with only 2 bonds, they act as surface rattlers. These rattlers further 
decrease the conductivity. Among d and e, e has lower conductivity because it has higher 
porosity. 

TABLE I. Thermal conductivity (TC), porosity, the ratio of number of atoms with 2-bonds to 
number of atoms with 4-bonds, (r2) and number of atoms with 3-bonds to number of atoms with 



4-bonds (r3) for the structures displayed in Fig 1101 

Struct. TC (W/mK) Porosity r2 r3 

a* 31.58 0.0127 - 0.0143 

b 23.55 0.0313 - 0.0265 

c 15.11 0.0313 - 0.0449 

d 6.54 0.0729 0.0159 0.0478 

e 5.75 0.0868 0.0162 0.0486 

f 8.47 0.0868 - 0.0823 



IV. CONCLUDING REMARKS 

In summary, we have developed a correct definition for energy moment for systems evolv- 
ing under classical dynamics. The thermal conductivity calculated from Einstein relation 
using the reformulated form of R yields the same result as the one calculated from HCACF, 
whereas the previous definition produces zero thermal conductivity in a non-diffusive solid 
system. We employed this method successfully for FCC argon and diamond silicon. Using 
-R and associated Einstein relationship simplifies the calculations compared to HCACF's. 
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FIG. 10. The atomistic arrangement of (a) 22-atom germanium clusters and (b)-(f) tubular va- 
cancies in silicon having different size and shapes. 



It overcomes the difficulties inherent in taking the analytical derivative of -R to obtain heat 
current. This also removes any ambiguity that arises from fitting the HCACF to a function 
in order to calculate thermal conductivity. Partition of R into potential and kinetic terms, 
moreover, demonstrates clearly the null contribution of potential energy term for any non- 
diffusive system. However, this is not case for systems with diffusion (directly contributing 
mainly as the convective heat flow). We also have demonstrated the utility of the method for 
bulk Si structures. With the studies performed on porous Si (cylindrical pores, and holes) 
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and Ge-nanodots embedded into Si matrix, we have explored specific features of R-based 
or HCACF based thermal conductivity calculations. In these systems, we observe increased 
noise on HCACF when compared to pure Si simulations, in contrast R-based calculation for 
these complex systems produced smoother behavior leading to well defined thermal conduc- 
tivity values. In our earlier studies, we have extensively taken advantage of this formulation 
in determining thermal conductivity of Si, Si-Ge nanostructures,— graphene nanoribbons 50 
and BN-nanostructures 51 . 
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